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1. Introduction 

In many applications one needs to integrate exact expressions involving several 
unspecified functions of multiple independent variables. The "homotopy operator" 
presented in this paper performs that task. Indeed, we give a formula based on 
a homotopy that integrates one-dimensional (1 D) expressions, or inverts a total 
divergence in two dimensions (2D) or three-dimensions (3D). The formulas are 
derived using a principle of the calculus of variations and are written in the language 
of standard calculus. 

The concept of homotopy, contributed to Poincare, is certainly familiar to differ- 
ential geometers and algebraic topologists P,'^. Briefly, a homotopy is a continuous 
map, T: X X [0, 1] — > Y between two functions u and uq (both are maps from a 
space X to a space Y), such that T(x, 0) = uo(x) and T(x, 1) = u(x). For example, 
T(x, A) = Au(x) -|- (1 — A)uo(x) = uo(x) -|- A(u(x) — uo(x)) is a homotopy. If uq = 
then r(x. A) = Au(x) expresses a scaling of u with a parameter A. The homotopy 
operator used in this paper is more elaborate for it integrates an expression with 
respect to an auxiliary variable A from to 1. 

Homotopy operators presumably first appeared in the inverse problem of the 
calculus of variations in works by Volterra They also appear in the proof of 
the converse of the Poincare lemma, which states that closed differential A:-forms 
are exact but only on domains with certain topological assumptions such as 
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a star-shaped domain [5]. Once the domain is appropriately restricted, the proof 
of the converse of Poincare's lemma requires the construction of a homotopy op- 
erator. More generally, the construction of suitable homotopy operators is the key 
to exactness proofs for many complexes such as the de Rham complex and the 
bi- variational complex [5]. 

One method to investigate the complete integrability 0| of a system of partial 
differential equations (PDEs) is to determine whether the system has infinitely 
many conservation laws. While developing a method for computing conservation 
laws of nonlinear PDEs in (1 + 1) dimensions 0, S], the authors needed a technique 
to symbolically integrate expressions involving unspecified functions. Mathemat- 
ica's Integrate function often failed to integrate such integrands, in particular, 
when transcendental functions were present. The computation of conservation laws 
of nonlinear PDEs in multiple space variables, requires a tool to invert total di- 
vergences 0, [13, [lH • The homotopy operator used in the proof of the exactness 
of the (bi-)variational complex [HI] can do the required integrations. The homotopy 
operator in [5i] has since been translated into the lang uag e of standard calculus and 
written into a more efficient algorithmic form [^, [gTllOl 11]. Although the homo- 
topy operator for one variable does not appear in work by Kruskal et al. [l^, it 
takes only one extra step to derive it. Taking that step, it also became clear how 
to generalize the formula to cover multiple variables. 

The calculus-based form given in [1, makes the homotopy operator easy to 
apply by researchers not familiar with differential forms. However, two issues must 
be addressed before the homotopy operator becomes a ready-to-use tool. The first 
issue relates to the inversion of a total divergence, Div^^, which does not have a 
unique answer. In analogy with standard integration with 'D~^, where the result 
is up to an arbitrary constant, Div~^ is defined up to curl terms. The homotopy 
operator will produce a particular choice for the curl terms [1, |l3], often creating 
very large vectors. Hence, we discuss an algorithm to remove the unwanted curl 
terms. Secondly, the homotopy operator fails to work on expressions involving 
terms of degree zero. A solution to this problem is also presented in this paper, 
thereby extending the applicability of the homotopy operator. 

With applications in mind far beyond the computation of conservation laws, 
the homotopy operator has its own Mathematica code, Homotopylntegrator .m 
(l3| . To our knowledge, Homotopylntegrator .m is currently the only full-fledged 
implementation of the homotopy operator in Mathematica. In collaboration with 
Hereman [o], Deconinck and Nivala [l3. [l5| have recently developed similar code in 
Maple. Although applicable to exact as well as non-exact expressions, their code is 
restricted to the 1 D case. Starting from version 9 of Maple, the homotopy operator 
is now part of the kernel of Maple to broaden the scope of its integrate func- 
tion. Anderson [l6| and Cheviakov 17, offer implementations of the homotopy 



operator in Maple, not as stand-alone tools but as a component of broader software 
packages. 

The multi-dimensional homotopy operator is an essential tool for the symbolic 
computation of fluxes of conservation laws of nonlinear PDEs involving multiple 
space variables. Most of such systems are not completely integrable for they only 
have a finite number of conservation laws [11]. As an application, we compute a 
conservation law of the (2 + 1)— dimensional Zakharov-Kuznetsov equation [l^ . 
which describes ion-acoustic solitons in magnetized plasmas. Recently, we devel- 
oped the Mathematica package ConservationLawsMD .m that automates the 
computation of polynomial conservation laws of nonlinear polynomial PDEs involv- 
ing multiple space variables and time. The multi-dimensional homotopy operator 
is a key tool in that package. 
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2. PreliminEiry Definitions 

Operations are carried out on differential functions, /(x, u'^^^(x)), where u^^^(x) 
denotes the dependent variable u = {u^, . . . , . . . , u^) and its partial derivatives 
(up to order M) with respect to independent variable x. Although we allow variable 
coefficients, the differential functions should not contain terms that are functions 
of X only. For simplicity of notation in the examples, we will denote n^, u^, u^, etc. 
by u, V, w, etc. 

We consider only ID, 2D, and 3D cases. Therefore, the independent space vari- 
able, X, will have at most three components. Thus, x represents x, {x,y), and 
{x,y,z) in ID, 2D, and 3D problems, respectively. Partial derivatives are denoted 
by Uk,xk2yksz = dfkl^d'ykTdX ■ For example, is written as U3x2y 

For simplicity, we often abbreviate /(x, u(^'*(x)) by /. The notation /[Au] means 
that in / one replaces u by Au, Ux by Aua,, and so on for all derivatives of u, 
where A is an auxiliary parameter. For example, if / = xu^Vx + UxV2x sin w then 
/[Au] = X^xu^Vx + X^UxV2xSinXw. 

Two additional operators are needed before the homotopy operator can be de- 
fined. The total derivative operator allows one to algorithmically compute deriva- 
tives of differential functions based on the chain rule of differentiation. 

Definition 2.1: With x = x, the ID total derivative operator T>x acting on 
/ = f{x,u^^\x)) is defined as 



j=l k=0 ^"'kx 

where is the order of / in component and M = maxj=i . . jyAfj'. The partial 
derivative, acts on any x that appears explicitly in /, but not on or any 
partial derivatives of . The total derivative operators in 2 D and 3 D are defined 
analogously. For example, the 3 D total derivative operator with respect to x is 

N Ml Mi Mi 

= ^ + E E E E "k+Da^ 



Qx ^ ^ ^ ^ {ki+l)xk2yk3Z r, j ' 
j=l fei=0A;2=0fe3=0 ""'kixk2yk3Z 

where M^, M2 and M3 are the orders of / for component with respect to x, y, 
and z, respectively 

Obviously, Vy and Vz can be defined analogously. Note that for the 3 D case, M 
is the maximum order for derivatives on , j = 1, . . . , N. 

The Euler operator, also known as the variational derivative, is one of the most 
important tools in the calculus of variations. 

Definition 2.2: The Euler operator for the 1 D case where / = f{x, u^^^\x)) is 



Mi 

^u^{x)f = y^(-^: 

A;=0 

df 



dui 



kx 



dui 



V: 



2x 



df_ 

^ dui 



+ ■■■ + {--Dx^'^ 



df 



3x 



dui 



Mix 



(1) 



j = 1, . . . ,N. The 2D and 3D Euler operators (variational derivatives) are defined 
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analogously. For example, in the 2 D case where x = (x, y) and / = /(x, u(^^)(x)), 
the Euler operator is defined as 

Ml Mi 

In this paper we use the Euler operators to test if differential functions are exact. 

Definition 2.3: Let / = /(x, u*^*^)(x)) be a differential function of order M. 
When X = X, / is called exact if there exists a differential function F{x, u(^^~^)(x)) 
such that / = D^F. When x = (x,y) or x = {x,y,z), f is exact if there exists a 
differential vector function F(x, u^*^"^) (x)) such that / = DivF. In the 2D case, 
F = and DivF = V^F^ + VyF^. The definition of Div in 3D is analogous. 

Theorem 2.4: A differential function f = /(x, u(*^)(x)) is exact if and only if 
'Cu(x)/ = 0- Here, is the vector (0,0, • • • ,0) which has N components matching 
the number of components of u. 



A proof for Theorem 12.41 is given in 11 1 



3. Integration by Parts Using tiie One-Dimensional Homotopy Operator 

In this section we show how to compute F = T)^"^ f = J f dx when / is exact. Direct 
integration by parts is cumbersome and prone to errors if / contains high-order 
derivatives of several dependent variables. The 1 D homotopy operator provides 
an alternate method for integrating exact differential functions (with one indepen- 
dent variable) and multiple dependent variables, often eliminating the need for 
integration by parts. 

Definition 3.1: Let x = x be the independent variable and / = f{x,u^^^\x)) 
be an exact differential function, i.e. there exist a function F such that F = T>~^ f. 
Thus, F is the integral of /. The 1 D homotopy operator is defined as 

where u = (u^, . . . , . . . , u^). The integrand, 1ui{x)f^ is defined as 

^.^(.)/ = E E-i(-^^)'-^^"'Mff' (4) 

k=l \i=0 / '^^kx 

where Ml is the order of / in dependent variable with respect to x. 

The usual homotopy operator with Aq = applies when ([3]) converges. However, 
due to a possible singularity at A = 0, ([3]) might diverge for Aq 0. This can occur 
with rational as well as irrational integrands. In such cases, one can take Aq — > oo 
or, alternatively, evaluate the indefinite integral and let A — > 1. As shown in [l5| . 
the homotopy operator with appropriately selected limits remains a valid tool for 
integrating non-polynomial functions. 

In the next section we prove that F = 'Hu(x)f- For now, we illustrate the appli- 
cation of the homotopy operator formulas ([3]) and ^ with an example. 
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Example 3.2 Let {u^,u'^) = {u,v) and take f{x,u^^\x)) = + 2xuux + UxV3x + 
U2xV2x — 3f^f2x- Applying the Euler operator ([T]) separately for u and v, with 
Ml = 2 and Mf = 3, one readily verifies that Cu[x)f = and L^(^x)f = 0- Thus, / 
is exact. To find F = J f dx, using first compute the integrands, 

df df 
'^u{x)f = ('"'-^)(^:^) + (""^'-^ ~ ^^''^^'dv~^ ^ u{2xu + V^x) + {uxl - uVx){v2x) 

= 2XV? + UxV2x, 



= + {vxl - vVx){^) + {v2xT - VxVx + vVt){^) 

= -v{QVxV2x) + {Vxl - vVx){u2x " 3v^) + {V2x1 - VxVx + v'DI){Ux) 
= UxV2x - ^vl, 

where I is the identity operator. Next, sum the integrands, replace n, Ux,Vx and V2x 
with Xu, Xux, Xvx, and Xv2x, respectively, divide by A, and integrate with respect 
to A. In detail, 

/■I 

F = T-(-u{x)f = {^u(x)f + 'Iv[x)f ) [Au] — 

= f {2Xxu^ + 2XuxV2x-^>?vl)dX 
Jo 



(A^xn^ + X'^UxV2x - >^^vl)\l 



= XV? + UxV2x - vl- 

Clearly, T^xF = f. Here Aq = was used since ([3]) converged for Aq — > 0. 

For polynomial differential expressions, ([3]) replaces integration by parts (in x) 
with a few differentiations following by a standard integration (of a polynomial) 
with respect to an auxiliary parameter A. The second example illustrates the inte- 
gration of a rational differential function where ([3]) diverges for Aq — > 0. 

Example 3.3 Let (u^jU^) = {u,v) and take f{x,u^^\x)) = {vux + uvx) / (uv)"^ . 
Using dH), Iu{x)f = ^ and Iy(x)f = Evaluation of gives 



F = Hu(x)f = I T3 — dX 



1 



Ao 



Obviously, Aq ^ would cause the integral to diverge. Instead, take Aq ^ 00 so 
that F = 'H^(^x)f — Alternatively, when a singularity occurs at A = 0, one 

could compute the indefinite integral and let A — s- 1. In either case, "DxF = f. 



4. The One-Dimensional Homotopy Operator 



The 1 D homotopy operator in Definition 13.11 is an expanded version of the one (in 



terms of higher-Euler operators) given in [9|, llJ, ll5|, ll8|] . Although they give the 
same result (ill . 21], the expanded homotopy operator formula is computationally 
more efficient for two reasons: the number of times that the total derivative is 
applied is significantly smaller and the combinatorial coefficients have been elimi- 
nated. 
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To prove that the homotopy operator ([3]) does indeed integrate an exact differ- 
ential expression, some extra definitions and theorems are needed. We first define a 
degree operator, A4 . Apphcation of to a differential function yields that function 
multiplied by its degree; hence its name. 

Definition 4.1: Let x = x be the independent variable. The degree operator A4 
acting on a differential function / = f (x , u^^'^\x)) is defined as 

N Ml 

j=l i=0 ^'^ix 

where / has order in with respect to x. 

The next example shows how the degree operator works. 

Example 4.2 Let u = {u^,u^) = {u,v) and /(x,u(^)(x)) = (u)^ (w2x)'^ (""5, 
where p, q, and r are nonzero rational numbers. Using ([5]), 



r 

X) 1 



Mf = U-^ {uf {v2xY {u5xY + U5x^— {uf {v2xT {u^xY +V2x^— {uf {v2xY {u^xY 
OU dU^x OV2x 

= pu {uY~^ {v2xY iu5xY + ^""50; {uY {v2xY (u^xY'^ + Q'V2x {uY {V2xy~^ {u^xY 

= ip + q + r) (uY {v2xY {u5xY ■ 

Note that the total degree of / has become a factor. Compare this with = 
x(nx"~"^) = nx" (n is a nonzero rational) in ID calculus. 

It is trivial to prove that is a linear operator. Less straightforward is finding 
the kernel (Ker) of A4. For what follows, we need the notion of homogeneity. 

Definition 4.3: A differential function / = /(x, u^^^^ (x)) is called homogeneous 
of degree p in u (and its derivatives) if /[Au] = A^/. 

Note the analogy with the definition of a homogeneous function of several vari- 
ables. For example for the two-variables case, /(x,y) is called homogeneous of 
degree p if /(Ax, Ay) = AP/(x, y). Examples are /(x, y) = ax^ + bxy + cy^ of degree 
2 and f{x,y) = Xj \J ax^ + hxy + cy^ of degree —1 (a, 6, and c are constants). 

Theorem 4.4: If f = j where k = k{x,u^^^\x)) and £ = £{x,u^^'^\x)) are 
homogeneous differential functions of the same degree (so, f is of degree 0) then 
Mf = 0. 

Proof : Suppose that k and i are homogeneous differential functions of degree p. 
Obviously, Mk = pk and M£ = pi. Hence, M{\) = ^i^^^-^i^^) = vihgihi = q. □ 

Conversely, one could ask "Which differential functions are in the kernel of the 
degree operator?" For simplicity, consider the case where f{u,v,Ux,Vx)- Any / E 
Ker M must satisfy Aif = 0, or explicitly, 

9f df df df ^ 
Utt + ^x^ ^'^^ + Vx^— = 0- (6) 

OU OUr OV OVr 



The jet variables u,v,Ux, and Vx are independent. Hence, (l6|) is a linear first-order 
PDE which can be solved with the method of the characteristic^]. Solving the 



^We are indebted to Mark Hickman for this argument and subsequent derivation. 
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first-order system of ODEs, 

du dux dv dvx 

dr or or ar 

where r parameterizes tlie characteristic curves, yields 

u = cie^, Ux = C2e^, v = cae^, Vx = c^e^ , (7) 

where the Cj are arbitrary constants. The first equation in d?]) determines e'^ = ^. 
Using it to ehminate r from the remaining equations of ([7|), yields 

C2 ^ C3 C4 

Ux = — u = Gi n, V = — u = L2 u, Vx = — u = C3 u, 

Cl Cl Cl 

where Ci,C2, and C3 are arbitrary constants. Thus, the general solution of ^ is 
/(Cl, C2, C3) = /(^, ^, ■^), where the functional form of / is arbitrary. 

Example 4.5 By construction, 

f = Bx (}^) = ^i--^--f) (8) 

is exact. Furthermore, Mf = by Theorem 14.41 Given f{u,Ux,v,Vx), not neces- 
sarily rational, one can a priori test whether or not / € KeiM. Indeed, applying 
the replacement rules, u ^ u, Ux — > fJ-u, v fiu, Vx — > fiu, to / should yield an 
expression which is independent of u, Ux,v, and Vx- Applying the replacement rules 
to / in dS]) gives -^r^, which only depends on fi. 

The argument above can be extended to any differential function / no matter 
which jet variables are present. The above "kernel test" can be performed on exact 
as well as non-exact differential functions. 

Next, the inverse of the degree operator will be derived, which requires the use 
of a homotopy [2I]. To avoid a potentially diverging integral at 0, the homotopy 
is defined on [Aq,!]. To begin with, consider only differential functions with one 
term, e.g., uux, xux, Ux/u,u^/ {^JvF+v^), or take / in Example 14.21 

Theorem 4.6 : Let f{x, u^^) (x)) be any single term differential function in 1 D 
and let 5(x,uW(x)) = Mf{x,u^^\x)), where Mf{x,u^^^\x)) / 0. Then, 

M-'g{x,n(''\x))= / ^[Au] ^. (9) 

Proof : If g{x, u*^^^)(x)) has order Mf in u^ with respect to x, then ^[Au] also has 
order in u^ . Furthermore, using the chain rule. 
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by the definition of Ad. Integrate botli sides of (jlOh with respect to A yields 

g[Xu]\l=M / 9[Xu]-, 
/•I J\ 

g{x,u^''Hx))-g[Xou] = M g[Xu]—. (11) 
To get dl]), ^[Aou] must be 0. This will affect the choice of Aq. Indeed, using ([5|), 

5[Aou] = EE^o<^^ = AoEE^^^- (12) 

j=ii=o 5(AoO j=ii=o 9{Xoui^) 

Depending on the form of /, there are two choices for Aq that make g[Aou] = 0. 

i. If ^[Aou] in ()12p is an expression in fractional form where Aq is a factor in 
the denominator, then let Aq ^ oo to get ^[Aou] = 0. The integral on the 
right-hand side of (llip will then converge for Aq ^ oo. 

ii. For all other cases, provided Aq does not drop out of (fT2l) after simplification, 
set Ao = to get (7[Aou] = 0. In the simplest case, when / (and therefore 
also g) is a homogeneous monomial of degree p, then ^[Au] has as factor A*^ 
and, consequently, g[Aou] = for Aq = 0. 

With g[Xou\ = 0, apply M'^ to both sides of ([U]) to get ([9]). □ 

Note that Theorem 14.61 excludes the case where M f {x , u'-^'^\x)) = 0. If 
/ G KerA^, then g{x,u^^^\x)) = and Theorem 14.61 becomes trivial. In view 
of Theorem 14.41 the homotopy operator cannot be applied to expressions of de- 
gree 0. In particular, the homotopy operator would incorrectly handle integrands 
involving ratios of polynomial or irrational differential functions of like degree. In 
Section [71 we provide a method to overcome this shortcoming of the homotopy 
operator. 

To show that ([9]) does indeed invert the A4 operator, in the next example we 
apply A4~^ to the result of Example 14.21 

Example 4.7 Let u = (n^,n^) = {u,v) and g{x,u^^'^^x)) = {p + q + r) 
(uf {v2xf {u^xY ■ Then, 

M-^g{x, u(*^)(x)) = [\p + q + r) {Xuf (A^^sx)^ (Aug.)^ ^ 



A 

= {p + q + r) {uf iv2x)' iu^xY C A^'+^+'-irf A 

Jo 

= (p + g r) [uf {v2xf {u5xY ( — : 

= {uf {v2xY {u5xY ■ (13) 

As with all polynomial integrands, we took Aq = as the lower limit of the homo- 
topy integral. Note that ()13p is identical to / given in Example 14. 2[ 

Remark 1 : Let f{x, u(*^)(j;)) = + H h /j H h /g, where each fi = 
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fi{x , u^'^-^^x)) is a single-term differential function (not necessarily a monomial). 
Theorem 14.61 holds for / provided that Mfi ^ for i = 1, . . . ,q. 

Next, we establish the commutation of the degree operator (and its inverse) with 
the total derivative operator. 

Lemma 4.8: MV^ = V^M and M'^V^ = V^M~'^ . 



Proof : A proof is given in ll[ . □ 



Now follows the key theorem for this Section which states that the 1 D homotopy 
operator does indeed integrate an exact differential function. 

Theorem 4.9 : Let f = f{x, u(^'^)(x)) be exact, i.e. V^F = f for some differen- 
tial function F{x,u'^^'^-^\x)). Then F = = 'Hu(x)f- 

Proof : The proof for the scalar case (u = u) in [21] can be generalized to the 
vector case. Indeed, first consider a fixed component of u and multiply Cu3(x)f 

■ df 

by to restore the degree. Subsequently, split off t{ — r. This yields 

ovJ) 

Mi Mi 

fc=0 '-^"'kx fc=l ^^kx 

■ df 

Next, integrate the last term by parts and split off u^. r. Continue this process 

dui 

df 

until j ■■ — is split off. Continuing with ([H]), the computations proceed as 

Mix 

follows: 



w 



df 



k=l kx k=2 kx 

Mi 



+ 4xY.(-^-' 



k-2 



Mix \ k=l kx 



k=2 ^^kx k=Mi "^"fcx^ 
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Mi (Ml -I Ml \ 

i=0 '^'^ix \ i=0 k=i+l ^'"'kx J 

2=0 ""-ix \k=l \j=0 / """fcx / 

Now, sum psp over all components to get 

N N Ml N (Ml /k-1 \ \ 

E = EE«i|f - E^. E E«L(-i'j'-"«' ^ 

j=l j=l i=0 ^^ii: j=l \k=l \j=0 / J 

Since / is exact, by Theorem 12.41 Ai^(t) f = for j = 1, . . . , A^, which implies that 
Ejii -Cu^Cx)/ = 0. Hence, 

/ N Ml /k-l \ ^, \ 

\j=l k=l \i=0 / '^'"fcx / 

Apply Jv[~^ to both sides and replace M~^Vx by VxM'^ using Lemma HTHl Thus, 
/ N Ml /fc_i \ \ 

/ = P. U-^ E E E -ic-^^o'^-^'^^^ ff • (17) 

\ j=l k=l \i=0 / "^^x J 

Apply V^^ to both sides of (fT7|) and use ([9]) to obtain 




The right hand side of ^ is identical to ^ with (jlj). □ 

The homotopy operator ([3|) has been coded in Mathematica syntax and is part 
of the package Homotopylntegrator .m (l3| . In extensive testing [lH, the ex- 
panded formula ([3|) out performed the implementation of the homotopy operator 
in [E, J^], dramatically reducing CPU time on complicated expressions. The code 
HomotopyOperator .m also integrates differential expressions that Mathematica 's 
Integrate function fails to integrate. For example. Integrate cannot integrate 
the (exact) expression 

f = UxV2xCOSU + V3xSmU - V4x- (19) 

The Homotopylntegrator .m code returns F = 'D~^f = V2x sinu — f3a;. It is easy to 
verify that "DxF = f. Needless to say, one does not need the homotopy operator to 
integrate simple expressions like (jl9p for it can easily be done with pen on paper. 
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The homotopy operator is a tool for integrating long and more complicated dif- 
ferential functions, in particular, those where Mathematica's Integrate function 
fails and integration by hand becomes intractable or prone to errors. 



5. Multi-Dimensional Homotopy Operators 

The 2D and 3D homotopy operators invert a divergence in 2D and 3D, respec- 
tively. Such inversions are considerably more difficult to do by hand. The 2 D and 
3D homotopy operators given in this paper were developed using the technique 
of proof in Theorem 14. 9[ Therefore, the homotopy operators presented below are 
different from the ones in 

0, E, Q, [ill, where they were defined in terms of 
higher-order Euler operators. Like in the 1 D case, the expanded versions of multi- 



dimensional homotopy operators require considerably less computations 

Definition 5.1 : Let f{x, y, u^^^^ (x, y)) be an exact differential function involving 
two independent variables x = {x,y). The 2D homotopy operator is a "vector" 
operator with two components, 

\'^u{x,y)J ' '^u{x,y)J J ' 

where 



(x) 

The x-integrand, XV |^^^^/, is given by 



(20) 



Ml Mi /fci-1 k2 \ Pif 

^1.)/ = E E E E^^^^<...(-^^)'^"^^"^-^.)'^"^0 d^'^''^ 

fci=l fc2=0 \ii=0 42=0 / '^^kixk2y 

with combinatorial coefficient B^^^ = B{ii,i2, ki, ^2) defined as 

(ii+i2\ /'fci+/c2-ii-«2-l\ 



B{h,i2,ki,k2, - ^^^^^^ 

r{y) 



Similarly, the y-integrand, '^^^^y-jf, is defined as 

Ml Ml / fci fc2-l \ 
fei=0fc2=l \ii=0 42=0 / ^%ia:A;2j/ 

where B^y^ = B{i2,ii,k2,ki). 

Definition 5.2 : Let /(x, u'-^''-* (x)) be a differential function of three independent 
variables where x = {x,y,z). The homotopy operator in 3D is a three-component 
(vector) operator, 

n_i{x) r n-ziy) f 1-1^^^ f 1 

' ^\i(x,y,z)J ' ' ^\i(x,y,z)J ' ' ^u{x,y,z)J I ' 
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where the x-component is given by 



The y- and z-components are defined analogously. The x-integrand is given by 



Ml Mi M^ fci-l k2 ki 
r{^) f ^ \- \- \- \- \- ^ r(x) J 



fci=l k2=0 ka=0 ii=0 12=0 43=0 

df 



'^'^kix k2y kiZ 

with combinatorial coefficient B^^^ = B{ii^i2,i^,ki^k2,k^) defined as 

Iii+i2+i3\ Ii2+H\ Iki+k2+k-i-ii-i2-ii-l\ Ik2+k-i-i2-i3\ 
I \ ki—ii — 1 I \ k2—i2 ) 



B{ii,i2,i3,ki,k2,k3 



^ki+k2+k3^ ^fca+fca^ 



As one might expect, the integrands I^^] ^ f and I^^} ^ f are defined analo- 

<=> f 1 o ui[x,y,z)-' u3(x,y,z)-' 

gously. Based on cyclic permutations, they have combinatorial coefficients B^^^ = 
B{i2,i3,h,k2,k3,ki) and B^^^ = i?(z3, ii, Z2, A^s, ^i, ^2), respectively. 

The homotopy with Aq = is used, except when singularities at A = occur. 
Theorem 5.3: Let f = f{x,u^^^\x)) be exact, i.e. f = DivF for some 
F(x,u(^^-i)(x)). Then, in the 2D case, F = Div^V = (^^(i j^)/' '^ufx j/)/) • 



ogously, ^n 3D, F = Div-^/ = (wlfi,,.)/' 

Proof : The very lengthy proof is similar to that of Theorem 14.91 It does not add 
much to the understanding of the theory. Nevertheless, details can be found in 11 



Appendix A] . □ 

The following 2D example demonstrates how (|20p works. 
Example 5.4 With u = {u,v) and x = let 

(3)r 



/(x, \V '{^)) = U Ux2y + 2uUxU2y + SUxVx COS W — AUyVxVxy — 2u2yV. 



X 



+V2y COS u + 3u2x sin v — UyVy sin u. (23) 

Using ([2T|) . compute 

(x) df df 1 df 

"(^■f) dux du2x 3 ^ ^ ^ ^ " 



x2y 



= 3n U2y + 3ux sinv, 
j(x) df I ^ \ 

I f \=V— (Vyl — VVy) 

^^(^'2^) dVx 2 ^ ^ dVxy 

= SUxVCOSV — 2UyVVxy — 2UyVxVy — 2U2yVVx- 
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Likewise, using ([22]) . compute 

I^H^ N = n-^^ + (uyl — uVy) — I — {2uxvT — UyT>x — UxT>y + 2uVx'Dy) 

"(^'2^) dUy ^ ^ ^' dU2y 3 ^ ^ ^ ^ dUx2y 

= —UVy sin u — 2uyv'^, 



I^f = + (y X - vVy) ^ + - ivxl - vVx) 7^ 

v{x,y) dy y y Qy ^ 2 ^ ^ dv 



y UV2y ^ UVxy 

■yV,j. -\- 2Uxyt 



2UyVV2x — 2Uyv'l + 2UxyVVx + Vy COS u 



Then, using ([20l) . compute 

^Iiis,)/ = {^u(x,y)f + '^^v(x,y)f ) I^"] " 

= / (3A^u^U2y + 3n^ sin All + 3Atia:f COS Af — 2A" 



UyVVxy 



-2X^UyVxVy — 2X^U2yVVx) dX 



U^U2y + ^Ux Sinv — "^UyVVxy — ^UyVxVy — \u2yVVx, 



dX 
~X 



'^u{x,y)f - i^u(x,y)f + ) [-^"1 

1 

[vy COS Xu — XuVy sin An — 4A UyVx + 2X UyVV2x + 2A UxyVVx) dX 

= ^UyVV2x + ^UxyVVx — ^Uyv'^ + Vy COS U . 

Thus, the homotopy operator gives the vector 

F = Div"^ f = ( "^""22/ + ^'^x si'^^ ~ l^yVVxy - luyVxVy - lu2yVVx \ ^24) 

y ^UyVV2x + ^UxyVVx — ^UyVx + Vy COS U J ' 

Obviously, there are infinitely many choices for F because the addition to F 
of a 2D "curl vector," K = (T>y9, —Dxff), where 9 is an arbitrary differential 
expression, will produce an identical divergence. Indeed, Div G = Div (F + K) = 
DivF + VxiTyyO) — T>y{'DxO) = DivF = /. The same happens in 3D where a curl 
vector, K = {T)yr] — T>zS,,'DzO — 'DxVi'^x£,~'^y^) involves three arbitrary differential 
functions 9,7] and ^. Often a concise result can be obtained for (j24p by removing 
undesirable curl terms. This is the subject of the next section. 

6. Removing Curl Terms 

The homotopy operator acting on an exact multi-dimensional differential function 
returns an F which includes a curl vector. Although the homotopy operator con- 
sistently leads to a particular choice for the curl vector, F can be of unmanageable 
size due to the presence of dozens, if not hundreds, of unwanted terms. Removing 
the curl terms makes F shorter. 

There are several methods for removing curl vectors. Doing so by hand requires 
educated guess work. Ideally, one could modify the integrand of the homotopy op- 
erator so that the least number of curl terms would be generated. We are currently 
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investigating this strategy. For now, we propose a simple, yet effective algorithm 
that uses only linear algebra. 

To remove the curl terms from (I24p . first attach undetermined coefficients, 
ki, . . . ,kp, to the terms of (I24p and call the new vector F. Doing so, 



kiu'^U2y + k2Ux sinv + k^UyVVxy + k^UyVxVy + k^U2yVVx 
kQUyVV2x + kjUxyVVx + ksUyVx + kgVy cos u 



(25) 



Next, compute 



Div F = 2kiuUxU2y + kiu^Ux2y + k2UxVx cos V + k2U2x sinv 

+ (^3 + ki + 2ks)UyVxVxy + {kz + k(i)uyVV2xy + (ks + kj) 
+ (A;4 + kQ)uyV2xVy + {ki + kj)uxyVxVy + (/cs + kQ)u2yVV2x 

+ {k5 + kj)Ux2yVVx + {k^ + k^)u2yv1 + kgV2yCOSU - kgUyVySiuU. (26) 

Since F in (I24p and F in ()25p should differ only by a curl vector, their divergences 
must be identical. Thus, DivF = DivF = / in ()23|) . Gathering like terms in (I23p 
and ()26p leads to the linear system 

ki = 1, k2= 3, /cs + A;4 + 2A;8 = -4, + fee = 0, 

fea + ^7 = 0, ki + ke= 0, ^4 + fey = 0, ^5 + fce = 0, (27) 

k5 + kr = 0, /C5 + fcg = -2, /cg = 1. 

To find the undetermined coefficients that will produce a F without curl terms, 
first solve for the ki that appear only in equations with non-zero right-hand sides. 
In ()27p . clearly ki = 1, k2 = 3, and kg = 1. The next variable to solve is kg since 
it appears in equations with non-zero right hand sides. Coefficients k^ through ki 
appear in equations with zero right-hand sides; these equations are solved last. 
Solving ([27ll in this order yields 

ki =1, k2= 3, A;3 = -kj, k^ = -kj, ^^g) 

k5 = -kr, ke = k-j, fcs = fey - 2, kg = 1, 



where k-j is arbitrary. Set fey = to eliminate as many terms as possible in (I25p . 
Substitute ^ with fey = into ([25]), to get 



-p ^ [ u^U2y + 3ux sinv , ^29) 

* Vy COS U — 2UyVx ' 



By construction, DivF = DivF = /, but (j29j) is free of curl terms. 

7. Extending the Applicability of the Homotopy Operator 

It follows from (jl6l) , that the 1 D homotopy operator in Definition 13.11 will not 
work when f{x,u^^\x)) (or for that matter any term of it) belongs to KerAi. 
Indeed, if TM/ = 0, in particular when / satisfies the conditions of Theorem 14.41 



then Ylf=i-^ui{x)f = C, where C is an arbitrary constant. Thus, the homotopy 
operator would return C, in fact for it ignores constants of integration. 

The next two examples (in 1 D) demonstrate a problem that occurs when some 
terms of / belong to Ker while others do not. 
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1 „,2^ 



Example 7.1 Let u = {u'^,u'^) = {u,v) and F = Then, 

U2xV - UxVx U2x UxVx 



f = VxF 



y2 y y2 



Obviously, / is of degree zero. Hence, A4f = by Theorem 14.41 Therefore, trying 
to integrate / with the homotopy operator will fail. Indeed, using @ for u{x) and 
v{x), respectively, yields 



r r , ( T T, \ Vx Vx Ux Ux 

Iu(x)f = Ut^ \- [uxl - uVx)t. = -u^ + u— H = — , 

^ ' OU-r OUOrr f V V 



T r df Ux Ux 
lv{x)J = = -V^ = 



Clearly, Iu{x)f + Iv{x)f = and ([3]) gives instead of F. 

Furthermore, if / is the sum of monomial, rational, or irrational terms (and frac- 
tions thereof) where any of these terms are in Ker A^, then the homotopy operator 
will return an incorrect result. Indeed, the homotopy operator would correctly in- 
tegrate the terms not in Ker M. , but annihilate the terms in Ker ^A . The following 
example highlights this situation. 

Example 7.2 Let u = {u^,u^) = {u,v) and F = Then, 
f = VxF 



U^Ux + U^Vx — 2uUxV + UVx — UxV 



{u — v)"^ 

U^Ux U^Vx UUxV UVx UxV , , 

+ 7 ^ - 2- ^ + ^ - , (30) 



(u — v)"^ {u — v)"^ {u — v)"^ {u — v)"^ {u — v)"^ 
is exact by construction. Applying the degree operator to /, term by term, gives 

U^Ux U^Vx „ UVUx 

= 7 + 7 ^ - ^7 12 +0 + 0- 

[u — w)^ [U — V)'^ [u — v)'^ 



Because the last two terms of (jSOp are annihilated by Al, the homotopy operator 
will incorrectly integrate f. One would obtain F = instead of F = " 

A method to overcome these shortcomings will now be proposed. A coordinate 
for the "jet" space where /(x, u(^)(x)) resides is given by 



u 



^^■^^(x) — {u^^ ,u\.,u\^, u\, u\^, u\y, u\^,u\y, . . . , U^jU^, . . . , U^N^ jyjNy M ^ ^ 



which has origin = (0, ... ,0). By shifting the coordinate away from 0, / can 
be taken out of Ker Ai. It suffices to shift one of the variables that appear in the 
denominator of the term in Kei A4. After integrating, the shift must be undone to 
put the integral at the origin. The following examples illustrate the procedure. 

Example 7.3 Returning to Example m| where f{x, u^^\x)) = "^"^"2""^" . Since v 
appears in the denominator, replace v with v — vq to get /o = _ Then 

use the homotopy operator to integrate /q. Using first compute 

Ut UtV 

-Lu(x)Jo = and -Lv{x)Jo = -7 yj- 

^ ' V — Vq ^ ' [V — Uo) 
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Next, using ([3]), compute 



^u{x)fo = / 7 i^uj — = - — dX 

Jo \v-vo [v-vqYJ X Jq (Xv-voy 



UxVo 



1 

Ux 



v{Xv - fo) V -Vo 
Finally, remove the shift by replacing v with v + vq (or simply set vq = 0) to get 

F = ^u(x)/ = — • 

Example 7.4 To integrate / from Example 17.21 the shift only needs to be applied 
to the last two terms of (I30p because these are in Kei A4. Therefore, we will apply 
the homotopy operator to 

U^Ux U^Vx UUxV {u — Uq)Vx UxV 



{u — v)"^ {u — vY {u — vY (u — uq — vY {u — uq — vY^ 

where u has been replaced hy u — uq in the last two terms. Note that it suffices 
to shift just one of the two variables in the denominator. Keeping the number of 
shifts to a minimum leads to a simpler integrand for the integration over A. The 
integrands for the homotopy operator are 

_ u^ju - 2v) uv _ u^v v{u - Up) 

^ ' [U — Uj^ [U — Uq — VY {U — V) [U — Uq — vY 

The first terms in each integrand have not been shifted. However, the second terms 
have been shifted. Applying ([3]), 

-xyW f _ /"V — 'm^__^\ 

2 1 



Xu ^ UqV 



U — V {u — v){Xu — Uq — Xv) 







u'^ I UqV V 



+ 



U — V \{u — v){u — Uq — v) {u — v) 



2 

U V 



+ 



U — V U — Uq — V 



Finally, setting uq = yields the desired result, F 



u +v 

U — V 



Similar difficulties might occur with functions in fractional form involving mul- 
tiple independent variables. Again, a shift of a coordinate in the denominator will 
solve the problem. However, the complexity of the integration over A increases 
rapidly when there are multiple independent variables (with shifts). Implementa- 
tion of the 2 D and 3 D homotopy operators (in Definitions 15.11 and 15. 2p requires 
the same amount of care. A discussion of the kernel of the homotopy operator, 
including an alternate strategy for dealing with the integration of homogeneous 



expressions can be found in 15l |. 
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8. Application: Computation of Conservation Laws of Nonlinear PDEs 

This Section covers the appUcation of the homotopy operator to the computa- 
tion of conservation laws of nonhnear PDEs. We developed a Mathematica pack- 
age, ConservationLawsMD .m [201], that automates these computations for nonlinear 
polynomial PDEs with a maximum of three space variables in addition to time. 

Definition 8.1: A conservation law for a given PDE with independent variables 
t and X is defined as 

Vtp + BivJ = 0, (31) 

where p = /9(t, x, u(*^)(t, x)) is the conserved density, J = J(t, x, u^^^ (t, x)) is the 
associated flux. T>t is the total derivative with respect to t and Div is the total 
divergence with respect to x. (j3ip is satisfied on solutions of the given PDE 0]. 

A conservation law is found by first computing the density, p, followed by the 
computation of the flux J. The latter will require the use of the homotopy operator. 

Following the approach by Hereman et al. [1, [^, , a candidate density is built 
as a linear combination (with undetermined coefficients) of differential terms which 
are invariant under the scaling symmetry of the given PDE. Once the form of p 
is determined, one computes T>tP and removes all time-derivatives using the PDE. 
By (fSTI) . Vtp must be a divergence. Thus, using Theorem 12.41 one requires that 

>C„.(x)(Ap) = 0, j = l,...,A^, 

where C^ii^^-^ is the Euler operator. This leads to a linear system for the unde- 
termined coefficients. Substituting its solution into the candidate for p gives the 
actual density. Finally, the flux J = — Div~^(P(p) is computed with the homotopy 
operator. 

To illustrate the method, we compute a conservation law for the (2 + 1)— dimen- 
sional Zakharov-Kuznetsov (ZK) equation [l^ which describes ion- acoustic solitons 
in magnetized plasmas. After re-scaling, the PDE takes the form 

ut + auux + /?(Au)x = 0, (32) 

where a and [5 are parameters and A = + is the Laplacian. 

8.1. Computation of the Density 

It is easy to verify that (j32p is invariant under the scaling (dilation) symmetry, 

(t, X, y, u) [X'h, X'^x, X'^y, X'^u). (33) 

This scaling symmetry can be computed as follows. Assume that ()32p scales uni- 
formly under 

{t, X, y, u) ^ (T, X, Y, U) = {X% X'x, X% A^, (34) 

with undetermined (rational) exponents a, b, c, and d. Applying the chain rule to 
(f32]l yields 

A'^-'^ Ut + aX^-"^"^ UUx + /^A^^"^ Usx + pX^^^^''^ Ux2Y 

= A"-'^([/t + aA''-^-'^ UUx + PX^^~^ Usx + /3A^+2'=-« Ux2y) = 0. (35) 
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One obtains ()32p in the new variables {T,X,Y, U), up to the common factor X°'~'^, 
ifb — d — a = 3b — a = b + 2c — a = 0. Setting b = —1, one gets a = —3, c = —1, and 
d = 2, which determines ([33 L ^ more algorithmic method for computing scaling 
symmetries can be found in a, [13] ■ 
The scaling symmetry carries over to conservation laws and therefore can be 



used to construct densities Indeed, since (f3T]) holds on solutions of ([32]) . both 
p and J must obey the scaling symmetry (133p of given PDE. For example, ()32p is 
a conservation law itself, with 

pW = n, = (ian2 + Pu2,, Pu,y), (36) 

expressing conservation of mass. It is straightforward to verify that 

= J(2) = (^n^ - p{ul - uD + 2/3m(u2x + ^X2y), -2/3u^%). (37) 



also are a density-flux pair of (j32p : the corresponding conservation law expresses 
conservation of momentum. Note that the densities in (j36p and (j37p scale with 
and A^, respectively. For brevity, we say that p^^^ has rank 2; p^'^^ has rank 4. The 
fluxes j(^) and J^^^ have ranks 4 and 6. As a whole, the conservation laws (|3ip 
with (/o^-^^ j(^)) and (/o^^^ J*^^-*) have ranks 5 and 7, respectively. 

We will construct a third density, p^^^ ^ which has rank 6, i.e. each monomial 
in p^^^ depends on u and its derivatives so that every term scales with A^. First, 
construct the list V = {u'^,u'^,u} containing powers of the dependent variable of 
rank 6 or less. Second, bring all terms in V with a rank less that 6 up to rank 6 by 
applying Vx and Vy. Considering all possible combinations of derivations yields 



Q = {u ,U^,UU2x,Uy,UU2y,UxUy, UU^y , U/^x , U^xy , U2x2y ,Ux3y, UAy } ■ 

Third, remove all terms from Q that are divergences (because they belong to the 
flux). For example, u^x = Div (lisa^ , 0) . This leaves 

Q = {u^,ul,UU2x,Uy, UU2y,UxUy,UUxy}- 

Fourth, find terms that are divergence-equivalent and remove all but the lowest 
order term. Two or more terms are divergence-equivalent when a linear combination 
of the terms is a divergence. For example, n^. and uu2x are divergence-equivalent 
since Ux + uu2x = Div {uux,0). Consequently, uu2x can be removed from Q. Doing 
so, Q = {u"^ ,Ux,Uy,UxUy}. If divergences and divergence-equivalent terms were not 
removed, they would lead to trivial and equivalent conservation laws, respectively. 

Now, form a candidate density by linearly combining the terms in Q with unde- 
termined coefficients Cj, 

P^^'' = ClU^ + C2Ux + CsUy + C^UxUy. (38) 



All, part, or none of the candidate density may be an actual density for (j32p . 
Indeed, the candidate density may be trivial or a linear combination of two or more 
independent densities. Computation of the undetermined coefficients will reveal the 
nature of the candidate density. 
With ([38]), compute 

Pfp(^) = SciU^Ut + 2c2UxUtx + 2c3UyUty + C4,{UyUtx + UxUty). (39) 
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Using ([32]) . replace ut with —{auux + Pusx + Pux2y)- Set E = —'Dtp^^\ to get 

E = SciU^iaUUx + Pu3x + PUx2y) + 2c2Ux{aUUx + Pu3x + PUx2y)x 

+ 2c3Uy{aUUx + Pusx + PUx2y)y + C4Uy{aUUx + /3u3a; + PUx2y)x 

+ C4Ux{aUUx + /?M3a: + PUx2y)y (40) 

Since (|10|) must be a divergence, use ([2]) and require 

^u{x,y)E = -2((3ci/? + C3a)u^U2j/ + 2(3ci^ + C3a)uyUxy + 3(3ci/? + C2a)uxU2x) 
+ 2ciaUxUxy + C4aUyU2x = 0. (41) 



It follows from (j4T]) that 

3ci/? + C2a = 0, 3ci/? + c^a = 0, ac4 = 0. (42) 

Before solving the system it is necessary to check for potential compatibility condi- 
tions on the parameters a and /?. This is done by setting each Cj = 1, one at a time, 
and algebraically eliminating the other undetermined coefficients. See, e.g., [ll, 22] 



for details about computing compatibility conditions. It turns out that there are 
no compatibility conditions for ()42p and the solution is 

C2 = -fci, C3 = -fci, C4 = 0, (43) 

where ci is arbitrary. Substitute (|13]) into ([38|) and set ci = 1 to get 

P^'^=n'-'-^{ul + nl). (44) 



8.2. Computation of the Flux 

After substitution of (gS]) and ci = 1 into (|iO]) . 

^ = 3n^(autij. + /3u3^ + /3'u^2j/) - ^Ux{auUx + Pu^x + I3ux2y)x 

- ^Uy{aUUx + I3u3x + (3Ux2y)y 

Since E = Div J^'^^, the flux J^'^^ can be computed with the 2 D homotopy operator 
which inverts divergences. Using (|2Up . the integrands are 

= + P {9u\u2x + lu2y) - 6ui3ul + uD) + ^ (6nL + 5uly + fuij, 

+ |m(u4j/ + ^2x21/) - Ux{l2U3x + 7Ux2y) - Uy{3u3y + 8U2xy) + \u2xU2y) , 

and 

^u{x,y)^ = ^Pu{uUxy - ^UxUy) - ^ (3-u(n3a3y + Ux3y) + Ux{lSU2xy + 3^3^) 
+ 5Uy{u3x + 3tij.2j/) - ^Uxy{u2x + U2y)) ■ 
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The 2 D homotopy operator formulas in (I20p yield 



uix,y) 1^ \ u{x,y) y L J ^ 

{SaX^u* + /3A^ (9m^(m2x + |w2y) - 6u(3u^ + u^)) 

- Uy{3U3y + 8U2xy) + ^U2xU2y)) dX 

= |an'^ + P {3u^{u2x + |i*2j/) - 2n(3u^ + Uy)) 

-Uy{\u:iy + ^U2xy) + \u2xU2y) , 

7^(f ,E = C (l^f ,e) [Xu] ^ 

= y (^3X^(3u{uUxy - ^UxUy) - ^X{3u{u3xy + Ux3y) +Ux{3U3y + 13U2xy) 

+ Uy{5U3x + 15na;2y) - \u2y{u2x + ^ixj/))) ^A 
= l3u{uUxy - AUxUy) - ^ {3u{u3xy + ^ajSy) + Uxil3U2xy + 3U3y) 
+ 5Uj;(m32: + 3Ux2y) - dUxy{u2x + U2y)) ■ 

The flux J = (n^^^^ y)^^^ulx y)^) a curl term, K = {VyO, -VxO), with 

,2„ , 



6* = 2Pu Uy + [3u{u2xy + USy) + 5{2UxUxy + 3UyU2y + U2xUy) 

After removing K with the technique in Section [6l the density-flux pair reads 

j(3) = (^^,4 ^ 3^^2^^^ _ g^^(^2 ^ ^2) ^ Mi(^2^ _ ^2^) _ ^[^xiu^x + t^x2|y) 

+ Uy{u2xy + ""Sy)), 3/3n^n^j^ + ^Ua;s/(li2z + ^^2^))- (45) 

There is one additional density-flux pair with explicit dependence on x and t, 



j(4) = (i(^^3 _ ^(^2 _ ^2) ^ 2/3w(n2. + U2y)) - lx{%u^ + /3n2.) + f n,, 

- 2[5(tUxUy + ^xuxy)) , (46) 

which expresses the conservation of center of mass. No other density- flux pairs 
than the four reported above could be found with ConservationLawsMD.m. 
Zakharov and Kuznetsov [3] found three "integrals of motion" which are identical 
to (p6]) . ([37l) . and (f46l) . Without mention of fluxes, Infeld [23] reports the densities 
in (I37p and (I45p and states that another constant of motion exists, but only for the 
ZK equation in 3D. Shivamoggi et al. 2j] claimed that there are only four conser- 
vation laws for the ZK equation, but gave different results. Shivamoggi computed 
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conservation laws for the potential ZK equation, vt + ^ + v^x + Vx2y = 0, derived 
from (j32p with a = /3 = 1 by setting u = Vx followed by an integration with respect 
to X. Doing so, he produced four nonlocal conservation laws for (I32p . The densities 
for the potential ZK equation reported in [2j] are correct but the fluxes for three 
of the four conservation laws are incorrect due to typographical errors [2^. 

In summary, the existence of only four conservation laws confirms that (|32p is 
not completely integrable. 
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